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Abstract. We consider an "elastic" version of the statistical mechanical monomer-dimer 
problem on the n-dimensional integer lattice. Our setting includes the classical "rigid" for- 
mulation as a special case and extends it by allowing each dimer to consist of particles at 
arbitrarily distant sites of the lattice, with the energy of interaction between the particles 
in a dimer depending on their relative position. We reduce the free energy of the elastic 
dimer-monomer (EDM) system per lattice site in the thermodynamic limit to the moment 
Lyapunov exponent (MLE) of a homogeneous Gaussian random field (GRF) whose mean 
value and covariance function are the Boltzmann factors associated with the monomer en- 
ergy and dimer potential. In particular, the classical monomer-dimer problem becomes 
related to the MLE of a moving average GRF. We outline an approach to recursive com- 
putation of the partition function for "Manhattan" EDM systems where the dimer potential 
is a weighted £i -distance and the auxiliary GRF is a Markov random field of Pickard type 
which behaves in space like autoregressive processes do in time. For one-dimensional Man- 
hattan EDM systems, we compute the MLE of the resulting Gaussian Markov chain as the 
largest eigenvalue of a compact transfer operator on a Hilbert space which is related to the 
annihilation and creation operators of the quantum harmonic oscillator and also recast it as 
the eigenvalue problem for a pantograph functional-differential equation. 

1 . Introduction. It is not uncommon that deterministic problems may reveal connections 
with stochastic counterparts, and it is particularly so if the original formulation inherently 
involves a pseudo-randomness production mechanism or rich combinatorics. When this 
happens, probability theoretic methods can provide an additional insight into the deter- 
ministic setting. Such are, for example, the problems concerned with chaotic dynamical 
systems under increasingly refined spatial discretization which models the finite computer 
arithmetic. Here, the effect of the round-off errors, magnified in the long run by positive 
Lyapunov exponents, leads to a subtle interplay between the ergodic behaviour of the un- 
derlying continuous system and the computer discretization. This can distort the invariant 
measure to the extent that, in place of its original absolute continuity, a persistent attract- 
ing centre emerges, thus destroying resemblance between the long-term behaviour of the 
discretised system and its continuous predecessor As reported in a series of papers by 
A.V.Pokrovskii, V.S.Kozyakin and their collaborators (see, for example, [39] and refer- 
ences therein), the attracting centre can be modelled through a finite Markov chain with 
an absorbing state. Moreover, a strikingly accurate numerical reproduction of statistics 
of discretised systems has been achieved with the subsidiary Markov chains obtained by 
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composing random maps with a common fixed point on the finite state space, sampled in- 
dependently according to a specially designed probability measure.^ The idea of random 
composite maps on a countable set can be extended to the continuum, for example, by con- 
sidering the compositions (pi o . . . o i^^. of jointly Gaussian real-valued random processes 
tfk{x) of a real parameter x. If these processes have differentiable sample paths and share 
a nonrandom fixed point a;,, then its stability under the maps t^^, as a random dynami- 
cal system, is determined (in the linear approximation) by the asymptotic behaviour of the 
products 

N 

■■= n (1) 

fe=l 

of jointly Gaussian random variables ^fe := as iV — > +00. Assuming that the 

processes (pk are all identically distributed, the standard ergodicity argument under a suit- 
able mixing condition (on the decay of correlations between the processes (pk with distant 
values of k) leads to the Lyapunov exponent 

lim ^=Eln|6|. (2) 

Here, the expectation E( ) involves only the common probability law of the processes ipk as 
if the latter were independent copies of a "mother" process, in which case (2) would follow 
immediately from the strong law of large numbers [43]. Similar products of appropriately 
rescaled random processes have attracted an active research interest due to their applica- 
biUty to multifractal modelling of turbulence cascades and other complex phenomena; see, 
for example, [1]. The problem of investigating the asymptotic behaviour of the products 
ipN from (1) becomes considerably harder if it is concerned with the moment Lyapunov 
exponent (MLE) which we define here as 

lim i^^, (3) 

provided this limit exists (otherwise, the upper limit or other limit points can be considered 
instead). In contrast to the conventional Lyapunov exponent (2), the MLE (3) is easy to 
compute only in the case of independent random variables in (1). If they are dependent, 
their joint probability law nontrivially enters the MLE even if correlations in the Gaussian 
random sequence ^1, ^2, Csi • • ■ decay fast enough to ensure mixing. Although the problem 
of computing the MLE as a modification to the usual Lyapunov exponent (2) may seem 
artificial in the context of random dynamics stability mentioned above, it has a remarkable 
connection with a statistical mechanical problem of computing the equilibrium free energy 
of a system of monomers and dimers on a lattice, and this connection is the main theme of 
the present paper. 

In the classical monomer-dimer problem [3, 11, 12, 13, 20, 25, 29], which is essentially 
solved in the planar case and remains unsolved in higher dimensions, the monomers are 
represented by singletons, and the dimers are two-element subsets of the lattice formed by 
nearest neighbours, with the energy of a dimer depending on its orientation. In particular, in 
a monomer-free "isotropic" setting on the plane, the partition function counts the number 
of "domino" tilings of a given finite subset of the planar lattice [3, p. 125]. It turns out 



^Note in passing that tlie randomizing effect of tlie spatial discretization was clarified in rigorous terms by the 
author of the present paper in [46] by showing that the round-off errors are independent and uniformly distributed 
(that is, form a non-Gaussian white noise-like sequence) with respect to a natural finitely additive probability 
measure (a frequency functional) on an algebra of quasiperiodic sets, provided the original dynamical system is 
nonresonant; see also [28, 47]. 
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that, for any dimension n, the number of such partitions of a finite subset A C Z" of the 
n-dimensional integer lattice coincides with the product moment 



MA(e) E n ex 



(4) 



of a specially constructed homogeneous Gaussian random field (GRF) := {(,x)xeZ"- 
Moreover, this GRF is of moving average type which is easy to simulate. This, in prin- 
ciple, provides an alternative probabilistic way (different from the probabilistic algorithm 
proposed in [27]) for approximate computation of the partition function in the monomer- 
dimer problem by modelling the auxiliary GRF and statistically estimating its product 
moments instead of the direct Monte-Carlo simulation of the underlying particle system 
which is complicated by the tilability issue. Due to this connection, the free energy of the 
monomer-dimer system per lattice site in the thermodynamic limit is proportional to the 
MLE of the auxiliary GRF, defined, similarly to (3), by 

ln|MA(C)| 



m(^) := lim 

A— i-oo 



#A 



(5) 



where is the number of elements in a finite set, and A — > oo is understood in the 
sense of van Hove [41]. This connection between the MLEs of homogeneous GRFs and 
the statistical mechanical partition functions remains valid for a wider class of monomer- 
dimer systems. 

In the present paper, we consider an "elastic" version of the monomer-dimer problem, 
which extends the original "rigid" formulation by allowing each dimer to consist of two 
particles at arbitrarily distant sites of the lattice, with the energy of interaction between 
the particles depending on their relative position. Thus, the geometric constraint that the 
dimers are formed only by nearest neighbours is removed, thereby introducing the long- 
range interactions; see Fig. I. In this more general setting, the partition function of the 




Figure 1 . Examples of configurations of the elastic dimer-monomer 
system in dimensions one, two and three (left to right). Monomers are 
depicted as "o"s. Dimer particles are represented by """s, whilst the 
bonds between them are shown as arcs in the one-dimensional case and 
as line segments in the higher dimensions. 

elastic dimer-monomer (EDM) system and its thermodynamic limit can be encoded in the 
product moments and the MLE of an auxiliary homogeneous GRF whose mean value and 
the covariance function are the Boltzmann factors associated with the monomer energy and 
the dimer potential, respectively. The structure of the covariance function singles out a 
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class of EDM systems where the dimer potential is a weighted ^i-norm (also known as 
"Manhattan norm") of the dimer vector A remarkable feature of the Manhattan dimer 
potential is that the corresponding auxiliary GRF is Markov and, moreover, belongs to the 
class of Pickard random fields (PRFs) [37, 38, 45]; see also, [6, 7]. PRFs were developed 
originally for the planar case and subsequently extended to three dimensions [22]. Unlike 
general Markov random fields, they are amenable to a unilateral single-pass simulation, for 
example, using the level-sets algorithm of [8]. The simplicity of the autoregressive-like 
modelling allows the Gaussian PRF (GPRF) to be utilized for the statistical estimation of 
the partition function of the Manhattan EDM system, similar to the role of moving average 
GRFs for the classical monomer-dimer problem. However, in addition to this utility for the 
alternative Monte-Carlo simulation, we take advantage of the specific Markov structure of 
the auxiliary GPRF to establish a recurrence relation for its conditional product moments 
which involves local transitions from one level set to another. Although the exposition is 
focused mainly on the two-dimensional case for better visualizability, a recurrence equa- 
tion for the conditional product moments of the GPRF can also be estabhshed in three 
dimensions. Note that similar recurrence relations are present in the statistical mechanical 
transfer matrix method for systems with short-range interactions [2, 3]. Since the interac- 
tions in the Manhattan EDM system are of long-range nature, the possibility of recursive 
representation of its partition function seems striking. 

For the one-dimensional Manhattan EDM system, where the dimer potential is pro- 
portional to the standard Euclidean distance, the auxiliary GRF becomes a homogeneous 
Gaussian Markov sequence satisfying an autoregressive equation. Due to the Markov prop- 
erty, the MLE is reduced to the largest eigenvalue of a compact transfer operator on an 
invariant cone in a Hilbert space which governs a recurrence equation for the conditional 
product moments. We obtain a matrix representation of the transfer operator in the ba- 
sis of Hermite polynomials and provide numerical results on the dependence of the MLE 
on two parameters of the univariate Manhattan EDM system. Although this eigenanalysis 
lies within the general framework of the transfer matrix method, we employ probability 
theoretic techniques and constructs such as measure change and conditional Gaussian dis- 
tributions. This allows the transfer operator to be linked with the annihilation and creation 
operators of the quantum harmonic oscillator [42, p. 91], and the eigenvalue problem to be 
recast into a pantograph functional-differential equation [5, 10, 26, 31, 44] which involves 
scaling of the independent variable. 

Note that the approach to the monomer-dimer problem through the product moments 
and MLEs of auxiliary GRFs, proposed in the present paper, is different from, and in fact, 
more general than, the classical methods based on Pfaffians of matrices (see [3] and refer- 
ences therein) and the links of the problem with the determinants of random matrices [30]. 
Our machinery differs also from the Gaussian integrals in Grassmann variables and their 
connections with the Pfaffians [18]. 

The paper is organised as follows. Section 2 describes the class of EDM systems being 
considered and formulates the problem of computing the equilibrium free energy in the 
thermodynamic limit. In Section 3, this problem is related to the product moments and 
the MLE of a homogeneous GRF. Section 4 shows that the auxiliary GRF for the classical 
rigid dimer problem is of moving average type. Section 5 relates the family of Manhattan 
EDM systems with GPRFs and focuses on the two-dimensional case. Section 7 obtains 
a recurrence relation for their conditional product moments associated with aggregated 
level sets of Section 6. Section 8 proceeds to one-dimensional Manhattan EDM systems. 
Section 9 derives the transfer operator which governs the recurrence equation for the condi- 
tional product moments of the auxiliary Gaussian Markov chain. Section 10 establishes an 
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inequality for the norm of this operator and a related upper bound for the MLE. Section 1 1 
considers the matrix representation of the transfer operator in the basis of Hermite polyno- 
mials and discusses its links with the ladder operators of the quantum harmonic oscillator. 
Section 1 2 relates the MLE for the univariate Manhattan EDM system with the eigenvalue 
problem for a pantograph functional-differential equation. Section 13 provides concluding 
remarks. 

2. Elastic dimer-monomer problem. Consider a system of identical particles residing 
at sites of a nonempty finite subset A of the n-dimensional integer lattice Z", with each 
site of A being occupied by exactly one particle. For some disjoint two-element subsets 
{x,y} C A, the particles at the lattice sites x and y are (chemically) bonded and form a 
dimer. The particles which do not participate in dimers are interpreted as monomers. 

The monomers and dimers are endowed with particular values of energy. The individual 
energies of monomers are all equal and their common value is denoted by V . The energy 
W{z) associated with a dimer with endpoints x ^ y is, that of the interaction between the 
particles in it and depends on their relative position z := x — y. Here, : Z" \ {0} — > 
K lJ{+°°} is ^ symmetric function which we will refer to as the dimer potential. 

The classical monomer-dimer formulation [12, 13, 25] allows dimers to consist of near- 
est neighbours only, so that all dimers have a fixed length and the energy of a dimer can only 
depend on its orientation. The "elastic" dimer-monomer (EDM) setting, outlined above, is 
more general. Indeed, it still leaves room for geometric constraints which can be imposed 
by putting W{z) := +oo for prohibited values of the dimer vector z. On the other hand, 
for admissible z, where W{z) is finite, the dimer potential W plays the role of a scoring 
function which favors dimers with smaller energy, though, in general, it does not forbid 
arbitrarily distant particles to form a bond. As in the classical rigid setting, the sites oc- 
cupied by monomers can also be interpreted as vacancies in which case the EDM system 
provides a lattice model for a gas of dimers. We will not, however, employ this alternative 
interpretation in what follows. 

A configuration lu of the EDM system is completely specified by a subset D of evenly 
many ~ 2r sites of A, and a partition {{xi,yi}, . . . , {xr, j/r}} of D into two-element 
subsets which represent dimers. The class of pair partitions of D is denoted by Vd, so that 
^f^Vo — (2r — 1)!!. With the complementary set A \ being occupied by monomers, the 
total energy of the particle system in the configuration lu is 

r 

U{uj):^{N-2r)V + Y,W{xk~yk), (6) 

k=l 

where N := #A. The class JIa {<^} of all possible configurations of the EDM system 
associated with the set A consists of 

LJV/2J , 

= E (^)(2r-l)!! = E((l + C)^) 

= 9fEe^(i+0|^^^ = aA'e^+AV2|^^^^,A^^^(_,) (7) 

partitions of A into two-element subsets and singletons. Here, [ J is the floor function, ^ 
is a standard normal (that is, zero mean unit variance Gaussian) random variable with the 
probability density function (PDF) 

i^{x) := e-^'/7v^, (8) 

and 

Hk{x) := (-l)'=e"'/29^e-"'/2 (9) 
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is the fcth modified Hermite polynomial. In (7), use is made of the even order moments 
E(C^'') — {2r — 1)!! and the moment-generating function Ee^^ = for the standard 

normal distribution in combination with the generating function of the Hermite polynomials 
[33, Theorem 2.5.6 on p. 7]: 

J2Hkix)y''/kl^e-y-y'^'. (10) 

fc^O 

By the variational postulate of Statistical Mechanics describing the canonical ensemble 
[21], the particle system, in contact with a heat bath at absolute temperature T > 0, acquires 
an equilibrium probability distribution over the configurational space 51a that delivers the 
minimum value 

Fa :=min ^ {U{u}) + TlnP{uj))P{uj) = -TlnZ^ (11) 

to the free energy functional over all probability mass functions (PMFs) P : Ha — > [0, 1] 
satisfying J^ugQa ~ ^' finiteness of the set f^A and the strict convexity of the 

functional being minimized in (1 1) over the simplex of PMFs, the equilibrium value Fa of 

A- 



the free energy is achieved at the unique Gibbs-Boltzmann PMF PA(t^) = e'^'^^'^^/Z, 



Here, 

/? := 1/T, (12) 
where the units of the temperature T are chosen so as to "absorb" the Boltzmann constant, 
and 

Za := J2 e"^"*"^ (13) 

is the canonical partition function [2 1 , 34] which encodes the equilibrium properties of the 
system. For example, the average total energy of the particle system is recovered from (13) 
as J^uenA ^('^)^a('^) ^ In^A- We will be concerned with the statistical mechan- 
ical problem of computing the partition function Za and the thermodynamic limit of the 
equilibrium free energy (11) per site of an unboundedly increasing fragment A of the lattice 
Z": 

li,, |A =_T lim (14) 

Although the convergence A — > oo is usually understood in the sense of van Hove [41], 
we will often assume, for simplicity, that A is a discrete hypercube {0, . . . , — 1}" (or 
another simply shaped subset of the lattice), where N is an integer parameter which tends 

to +00. 

3. Auxiliary Gaussian random field. Suppose the temperature T is fixed and the dimer 
potential W is extended to the origin so as to satisfy the condition 14^(0) 7^ 00 and the 
inequality 

«eZ"\{o} 

To make this possible, we assume that the sum on the right-hand side of (15) is finite. Then 
the function C : Z" — > M_|_, defined as the Boltzmann factor associated with the dimer 
potential W by 

C{z):^e-^^'^'\ (16) 
with the shorthand notation (12), is positive definite [15]. That is, for any positive integer m 
and any points 21, . . . , in Z", the matrix {C{zj — is positive semi -definite. 

The fact that (15) implies the positive definiteness of C in (16) is established by a diagonal 
dominance argument [19, p. 349]. Hence, by Bochner's theorem, C is the covariance 
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function of a real-valued homogeneous random field on Z" which can be chosen to be 
Gaussian (without additional assumptions on C). Therefore, the fulfillment of (15) ensures 
the existence of a homogeneous Gaussian random field (GRF) ^ {£,x)x£Z" satisfying 

E^. = e-''^, cov(^,, Cy) = C{x - y) (17) 

for all x,y G Z". Here, cov(-, •) is the covariance of random variables, V is the monomer 
energy, and C is given by (16). In what follows, such a random field ^ will be referred to as 
the auxiliary GRF. Its significance for the EDM system, described in Section 2, is clarified 
by the theorem below. 

Theorem 3.1. Suppose the dimer potential W satisfies (15). Then for any finite set A C 
Z", the partition function (13) coincides with the product moment (4) of the auxiliary GRF 
^ over A: 

Za = Ma(0. (18) 
Proof. The product moment of the auxiliary GRF ^ over the set A can be computed as 

LJV/2J 

Ma(0 = l^''-*''Mn{v) = E E (19) 

DCA r=0 DcA:#D=2r 

where N := ^A, and rj := {rix)xeZ" is a GRF obtained by centering ^ as rj^ := — fJ', 
with 

■■= e-^^, (20) 
in conformance with (17). By the Isserlis-Wick theorem [23], which relates mixed mo- 
ments of evenly many zero mean Gaussian random variables with their cross-covariances,^ 
it follows that 

r r 

^d{v) = e n = E n - yk). (21) 

Vo k=l Vd k=l 

Here, the sum extends over the class Vd of partitions {{xi,yi}, . . . , {xr, yr}} of the set 
D into two-element subsets, with r :— #Z?/2. Now, upon substitution of (16) into (21), 
(19) takes the form 

[N/2i I r \ 

MA(e)= E e-f^-^"-)^^ E exp U/^Em^/c-y.) 

r=0 _DcA:#D=2r V k=l ) 

E E^^P ((^-#^)^+ E ^(^^-2^'=)) I • ^22) 

DcA:#Z?isovon Vo \ \ k=l J J 

In view of (6) and the definition of the configurational space JIa of the EDM system from 
Section 2, the right-hand side of (22) coincides with the partition function Za in (13), thus 
proving (18). □ 

In view of Theorem 3.1, the rightmost limit in ( 1 4) coincides with the moment Ly apunov 
exponent (MLE) of the auxiliary GRF ^: 

A-i-CO ^A A-i-OO :^A 

which is identical to (5) since the product moments of the GRF being considered are all 
nonnegative. This limit is completely specified by the mean ^, given by (20), and the 



^and is used as an algebraic moment closure in tlie statistical tlieory of turbulence [14, pp. 44^5] and quantum 
field theory (see also [24, Theorem 1.28 on pp. 1 1-12]) 
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spectral density function S : [— tt, tt]" defined as the Fourier transform of the 

covariance function (16) by 

S{X) J2 Ciz)e-"^"' = ^ e'''^^^^-^^'^ (24) 

Here, m"""?; is the inner product in M", vectors are organised as columns unless indicated 
otherwise, and (•)^ denotes the transpose. Computing the MLE m(^), with Ma(0 in 
(23) being replaced with |Ma(C)|i is of considerable interest in its own right for a general 
homogeneous GRF ^. However, the salient feature of this problem in the context of the 
EDM system is that, in view of (16), both the mean value and the covariance function of 
the auxiliary GRF in (17) are nonnegative. The nonnegativeness of C implies that the spec- 
tral density S in (24) is not only symmetric and nonnegative, but is also positive definite. 
Hence, S is the covariance function of yet another homogeneous GRF, this time, on the n- 
dimensional torus [— tt, tt)", whose probability distributions are invariant under the group 
of entrywise additions modulo 27r. 

Remark 1. The fact that the choice of 1F(0), which affects C(0) in (16), has no influence 
on the product moments of the auxiliary GRF ^ is explained by their invariance under 
adding a ^-independent "white noise" GRF ( {Cx)xeZ'^ whose values are independent 
Gaussian random variables with zero mean and variance j'^. The resulting GRF ^ + :— 
+ Cx)xei,^ has spectral density S" + 7^ and its product moment over a finite set A C Z" 
is computed as 

MA(e + = EMA(e + c Ua) = E n (e- + EC.) - MA(e). (25) 

xeA 

Here, ^a := {£,x)x€A is the restriction of ^ to the set A, and, in extending (4), use is made 
of the conditional product moment Ma(?/ | ©) := E (HajGA Vx | ©) of a random field 
f] := {r]x)x£Z'^ over the set A with respect to a ci-subalgebra & (or a conditioning random 
element which generates &). In (25), we have also used the tower property of iterated 
conditional expectations [43], and the assumption that C consists of independent zero mean 
random variables, independent of f . Therefore, under the transformation jj. ^ pfi and 
S I— >■ p^S + j'^, which corresponds to + C with a constant factor p > 0, the MLE 

m(^) in (23), as a functional of p and the specti'al density S from (24), is transformed as 
m > m + In p. 

4. The classical monomer-dimer problem and moving average Gaussian random fields. 

Consider the structure of the auxiliary GRF for the classical monomer-dimer problem 
[12, 13, 25] as a particular case of EDM systems. Suppose the dimer potential W on 
the set Z" \ {0} is given by 

Wiz):={ ^ fr^i'' > (26) 

^ ' \ +00 if |z| > 1 ' ^ ^ 

where ak is the energy ascribed to dimers which are parallel to the fcth coordinate axis, and 

efe (0^_^,1, 0^_^)T (27) 

k — l zeros n—k zeros 

is the fcth basis vector in M". Note that (26) forbids dimers {x, y} with |x — y| > 1 and 
describes the dependence of the energy of a dimer of Euclidean length 1 on its orientation. 
In particular, ifV:= +oo (so that monomers are energetically forbidden) and «! = ...= 
a„ = 0, then the partition function Z\ in (13) in this monomer- free isotropic case counts 
the number of "domino" tilings [3, p. 125] of the set A, that is, partitions of A into pairs 
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of nearest neighbours, assuming that ^^A is even. In general, the dimer orientations with 
smaller values of ak are energetically favoured. The resulting series on the right-hand 
side of (15) consists of finitely many nonzero terms, and the auxiliary GRF ^ for this 
classical "rigid" dimer system can be explicitly constructed as follows. Using (27), for 
every k = 1, . . . , n, we denote by 

^fc:=Z" + efe/2 (28) 

the set of midpoints of the edges between neighbouring sites in Z" parallel to the fcth 
coordinate axis. For any point x £ Z", the nearest sites of Ek are x ± ek/2, both at the 
Euclidean distance 1/2 from x. Accordingly, those midpoints from the combined set 

n 

E:=\jEk, (29) 

fc=i 

which are nearest to a given x e Z", are x ± 6^/2, with k = 1, . . . ,n. We will now 
associate independent standard normal random variables rjy with the sites y E E, and 
ascribe, to every x :— {xk)i^k^n G Z", a weighted sum 

n 

:= M + X! '/Pk{'nx~ek/2 + Vx+ek/2), (30) 

fc=l 

over the 2n sites of the set (29), nearest to x, where pi , . . . , p„ are nonnegative parameters 
given by 

Pk:=e-">'^. (31) 
The resulting GRF ^ :— {£,x)x£i," is homogeneous and has mean /i and covariance function 

C{z) := cov(^o,Sz) = < Pk if 2 = ±efc , (32) 

otherwise 

which, in view of (31), is related with the dimer potential (26) by (16), with W being 
extended to the origin as W{0) = —T\n{2jy^^i Pk)- Therefore, (30) indeed produces 
an auxiliary GRF ^ for the classical dimer system. It is a moving average random field 
which is easy to simulate. The simulation boils down to generating the standard normal 
random variables at sites of the set E given by (28)-(29). This, in principle, provides an 
alternative probabilistic way for approximate computation of the partition function Z\ for 
the monomer-dimer problem over a finite set A C Z" by modelling the auxiliary GRF ^ and 
statistically estimating its product moment Ma(^) instead of the Monte-Carlo simulation 
of the underlying particle system which is more complicated. 

In view of (32), the covariance matrix cov(^a) of the restriction of ^ to a finite set 
A C Z", with being considered as a Gaussian random vector of dimension #A, has a 
(2n+l)-diagonal structure. A similar sparsity structure is known for the inverse covariance 
matrices (that is, precision matrices) of finite fragments of Markov GRFs [40]. The Markov 
GRFs describe the equilibrium states of systems of linear harmonic oscillators with nearest 
neighbour interaction and are, in general, hard to simulate in contrast to the moving average 
random fields above; see also [4, 16]. An important subclass of Markov GRFs which is 
amenable to efficient unilateral simulation in a single pass through the simulation domain 
is provided by Pickard random fields (PRFs) [6, 7, 22, 37, 38, 45] which behave in space 
like autoregressive processes do in time. 
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5. Manhattan dimers and Gaussian Pickard random fields. Suppose the dimer poten- 
tial is given by a weighted £i-norm (also called "Manhattan norm") of the dimer vector 

n 

:=^afc|zfe|. (33) 

k=l 

Here, ai, . . . , a„ are positive parameters which quantify the projections of the force of 
attraction between the particles in a dimer onto the coordinate axes. Unlike ideal springs, 
the hypothetical Manhattan dimer does not develop a linearly increasing force when being 
stretched, albeit the attraction between the particles in it persists. The covariance function 
(16), associated with the Manhattan dimer potential (33), is 

n 

fc=i 

where < pi, . . . . p„ < 1 are computed according to (31). The resulting auxiliary GRF 
^ is a Gaussian PRF (GPRF) in dimensions two [45] and three [22] and lends itself to the 
unilateral simulation mentioned above. Such simulation of the GPRF ^ is made possible by 
the following enhancement of the spatial Markov property for any x :— {xk)i^k^n G Z": 

ICa-1 = K. I^Aol. (35) 

Here, fr] | is a shorthand notation * for the conditional probability distribution of ij given 
and 

A- := Z"\A+ A+ Z" n x [xk, +oo), A^ f x {xk - 1, Xk}) \{x}. (36) 

The sets A^ and A+ are the "past" and "future" of the lattice with respect to the site x in 
the spatial sense [7]. Accordingly, the set A", which consists of 2" — 1 points from A^, 
plays the role of the nearest past for x; see Fig. 2. In the one-dimensional case n = 1, 
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Figure 2. A fragment of the two-dimensional integer lattice Z'^. The 
sites of the spatial past A^ and the spatial future A+ for a given site 
X e 7?, defined by (36), are represented by "o"s and """s, respectively. 
The sites of the nearest past A|^ are marked by "-k^'s. 

which we will study in Sections 8-12, the relation (35) between the conditional probability 



•^We use I?; | (^| instead of the more standard notation P^|^ in order to mitigate tlie burden of multitiered 
subscripts. 
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laws reduces to the standard Markov property, thus allowing ^ (which becomes a Gaussian 
random sequence) to be generated by an autoregressive equation. 

A remarkable feature of the Markov structure of the GPRF in the multivariate case 
n > 1 is that its restrictions to the level sets Lq, Li, L2, ■ ■ ■, described below, form 
a Markov chain of order n. This property is employed by the level-sets algorithm of [8] 
for the unilateral simulation of such random fields on the nonnegative orthant Z" of the 
lattice, where Z_[_ is the set of nonnegative integers. For example, in the two-dimensional 
case n — 2 elucidated by Fig. 2, the level sets are given by 

:= {{j, k)&Zl: j + k^e}, £ ^ 0, (37) 

and the restrictions of the GPRF, which are organised as Gaussian random vectors 
(with (_Lg having dimension £ + 1), form a Markov chain of order two. More precisely, 
for any £ > 0, the conditional probability distribution ICl^+i | £.Lo, ■ ■ ■ ,£,Lel coincides 
with I £,Le-iJ ^LeJ- Furthermore, the values of ^ at sites of the next level set i^+i 

are conditionally independent given ^Lf_i and ^l^, with their conditional joint probability 
distribution 

(j,k)eL!+i ok 

being the product of £ + 2 conditional Gaussian distributions. Here, at site (j, fc) is 
conditioned on the values of ^ in the nearest past A°j, which is restricted to the nonnegative 
quadrant of the lattice as 

r {0--l,fc-l), (i,fc-l), 0--l,fc)} fori>0, fc>0 
AO, := K% n ^+ ^ S {(j - 1, 0)} for j > 0, fc = 

[ {(0,fc- 1)} for j = 0, fc > 

(39) 

in conformance with (36) for the two-dimensional case. In terms of the centered random 
field (T := {(Jx)x&'^, defined by 

(Jx ■■= L - M, (40) 

with /i associated with the monomer energy by (20) as before, the conditional Gaussian 
distributions on the right-hand side of (38) are given by 

{J\f{aaj^i^k-i + b(Jj.k-i + ccrj-i,fe, d^) for j > 0, fc > 
AA(pia,_i,o, 7?) forj >0, fc = , (41) 

A/'(p2CTo.fc^i, 7I) for i = 0, fc > 

where pi, p2 are related to the attraction force parameters ai, a2 of the Manhattan dimer 
potential (33) by (31), and the conditional variances 7J, 7I of the last two distributions are 
computed as 

71 := ^1-pI 72 := ^ I - pI (42) 

This follows from the well-known results on conditional distributions for jointly Gaussian 
random variables [43]. The coefficients a, 6, c and the conditional variance (P of the 
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Gaussian distribution for j, fc > in (41) are computed as 



[a, b, c] 



= [PlP2, P2, Pi] 



1 

Pi 
P2 



Pi 
1 

PlP2 



P2 
PlP2 
1 



= [-PlP2, P2, Pi] 



(43) 



var((Tjfc) - [a, 6, c] cow{ajk,^]k)^ 



1 - [-PlP2, P2, Pi] 



PlP2 
P2 
Pi 



(1-P?)(1-P2) 



2 2 

7i72, 



(44) 



with var() the variance of a random variable. Here, the conditioning restriction ctto of the 

centered random field (40) to the set (39) is organized as the three-dimensional row-vector 

■ujjk '■= [(Tj_i_fc_i, (Tj.fc_i, CTj-i^fc] whose covariance (3 x 3)-matrix cov(ci7jfc) and the 
row-vector cav{crjkT^jk) of cross-covariances with ajk are appropriate submatrices of 
the covariance matrix 





^3-l,k-l 
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Pi 
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Pi 
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Pi 
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_ £,jk _ 
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PlP2 


P2 


Pi 


1 



obtained by evaluating (34) for the square cell {j 
also [38, 45]. 



1, j} X {fc — 1, fc} of the lattice; see 



6. Aggregation of level sets. We will now employ the technique of reducing the order of 
a Markov chain by augmenting its state space. For any £ > 0, consider the restriction of 
the centered GPRF cr, defined by (40), to the union 



Le-i [J Li 



(45) 



of two adjacent level sets from (37) which consists of 2£+ 1 lattice sites. Also, let Oo '■= Lq 
be the singleton formed by the origin of . The resulting sequence of random vectors ctq^ , 
labeled by ^ ^ 0, is a Markov chain of order one. In view of the reflection symmetry of 
the covariance function (34) with respect to the coordinate axes, the restrictions C-ef of 
a to the sets —6^, which are centrally symmetric to the aggregated level sets (45) about 
the origin, also form a Markov chain. Furthermore, the Markov property remains valid 
upon removal of the "endpoints" (0, —£) and (—^,0) from the set —Qg. More precisely, 
the restrictions (Tr^ of the GPRF a to the sets Fq, Fi, . . . defined by 

r^v (-eAr+i)\{(0,-A^-l), (-7V-1,0)} 

= {(0, -N), (0, -N - 1), (-1, -A^ - 1), . . . , 

{-N + 1, -1), i-N, -1), (-A^, 0)}, (46) 

also form a Markov chain. We label the values of a at sites of the set F^r by integers to 
2N as shown in Fig. 3, which allows (Tr„ to be identified with a {2N + 1) -dimensional 
Gaussian random vector 



N,k)0<k<2N 



■— {o'O.-N, CTO-N-l, <^-l-N-l, 
<^-N+l-l, <^-N,-l, 0--N,o}- 



(47) 



Due to this numbering, Sat is obtained by sampling the GPRF a along a monotone path 
in the planar lattice. Hence, by the results of [38, 45], for every given A ^ 0, the entries 



MOMENT LYAPUNOV EXPONENTS OF GAUSSIAN RANDOM FIELDS 



13 



-N 



O O O ( > 



o o 
-® o 



-N 



Figure 3 . The subset Aat of the nonpositive quadrant of the planar 
integer lattice, defined by (52), is represented by "o"s for iV — 4. The 
sites of the set F jv from (46) are marked by "★"s and labeled by integers 
to 2N in a "zigzag" manner (along the dashed line) starting from the 
point (0, -A^). 



'^Nfl, ■ ■ ■ ,^N.2N of the random vector S n form a nonhomogeneous Markov chain. More- 
over, by evaluating the covariance function (34) for the zigzag shaped set T^, it follows 
that Eat is an alternating autoregressive sequence [36] governed by the recurrence equation 

„ _ J pi^N.j +JiWn,j for j even 

^""''^'^X P2^NJ+72WN,, for J Odd ' ^^^^ 

where the initial element Sat^q and the innovations wj^^, wn,!, wn.2, • • ■ are independent 
standard normal random variables, and use is made of (42). For any iV > 0, the conditional 
Gaussian distribution of Y^n-i given S^v is degenerate since these random vectors have 
common entries: 

YN,2k = ^N-iak-u l^k<N. (49) 
The other entries of Sjv_i (namely, SAr_i,2fc for ^ fc < N) are conditionally indepen- 
dent given Etv- More precisely, 

f{YN-1.2k)o^k<N I SatI — X |I]Af_i.2fc | 'YN.2k, Sjv,2fc+1, YN,2k+2j, (50) 

A;=0 

where each of the conditional distributions on the right-hand side is described by the first 
of the Gaussian distributions from (41): 

|SjV-l,2A: I YN,2k, Sjv,2fc+1, Sjv^2fe+2l 

= N{p2YNak - PlP2YN,2k+l + PlYN,2k+2, ■ (51) 

Remark 2. Notice the local nature of the conditional distribution iS^r-i | SatJ. Indeed, 
(49)-(50) imply that IEat-ij | Sa'1 = | Satj, Satj+i, SAf,i+2l involves at 

most three entries of the random vector E for every ^ j ^ 27V — 2. This property is in- 
herited from the localness of the conditional probability measures describing the transitions 
of PRFs between the level sets. 

7. Recurrence relation for conditional product moments. Assuming > and using 
the sets Fq, . . . , Fat from (46), we will now consider the two-dimensional Manhattan EDM 
system on the set 

:= U r,, (52) 
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which consists of #Ajv = N{N + 5)/2 + 1 lattice sites; see Fig. 3. The product moment 
of the auxihary GPRF ^ over A^r is representable as 

Ma„(0 = Ei?jv(Sw). (53) 
Here, the function Rn : — > M is defined as the conditional product moment of ^ 
over the set A^r with respect to the random vector Ejv from (47): 

i?Ar(I]jv) :=MA„(e I Sat). (54) 

In particular, since the set Aq is a singleton which consists of the origin of Z^, the definition 
(54) yields 

i?o(s) =E(Coo I croo = s) = M + (55) 

Theorem 7.1. For any N > 0, the conditional product moment function Rn of the aux- 
iliary GPRF ^ for the two-dimensional Manhattan EDM system satisfies the recurrence 
equation 

Rn{^n) = E(_Rjv_i(S]Ar_i) I Ejv) 

X{fi + Sw^o)(A^ + ^N,2n) H (M + ^N,2k+l), (56) 

k=0 

with the initial condition (55). Here, the conditional expectation is taken over the condi- 
tional probability distribution |SAr_i I StvI described by (49)-( 51). 

Proof. By partitioning the set from (52) into two disjoint subsets AAr_i and \ 
Ajv-i — Tjv \ Ajv-i, with Fjv given by (46) (see also Fig. 3), it follows that the function 
(54) can be factorized as 



RNi^N) = E n u 

= MA„_,(e|SAr) n (57) 

a:erjv\Ajv-i 



Here, 



Y\_ Ca: = + SjV,o)(Ai + SAr^2Ar) J]^ (M + SiV,2fc+l) (58) 
2;erjv\Aiv-i fe=o 
is a polynomial function of the appropriate entries of the random vector S^y in (47). Fur- 
thermore, 

Ma„_,(^|Sjv) = E(MA„_,(e|Sw_i,EAr)|I]Ar) 

= E(i?jv-i(Siv-i) I Sw), (59) 
where, in addition to the tower property of iterated conditional expectations, we have used 
the relation |(^ An- 1 I ^A'-ii^Jvl — Kajv-i I SAr_i| which follows from the measur- 
ability of £,An-i with respect to Eq, . . . , Sat^i and the Markov property of the sequence 
Soi • • • , Sjv. The recurrence equation (56) is now obtained by substitution of (58) and (59) 
into (57). □ 

By induction, (55) and (56) imply that the function i?Ar is a polynomial (in 2N + 1 
scalar variables) of degree ^A^r- Therefore, calculating the partition function of the two- 
dimensional Manhattan EDM system for the set A^ through the product moment of the 
auxiliary GPRF ^ according to (53) consists in evaluating the expectation of the recursively 
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computed polynomial i?jv on the segment Ejy of the alternating autoregressive sequence 
governed by (48). For example, by recalling the conditional distribution (41)-(44), it fol- 
lows that 

-Ri(so, si, S2) = E(_Ro(croo) I cTo-i = So, = Si, o--i,o = S2) 

+ Sq){^ + Si){^ + S2) 
= (M + P2S0 - P1P2S1 + PlS2){^i + so)(m + + S2) 

is a quartic polynomial. Further iterations of the recurrence relation (56) involve the con- 
ditional expectations of higher order polynomials applied to conditionally independent 
Gaussian random variables; cf (50), (51). Although this recurrence is not easy to im- 
plement, its significance is that the above discussed localness of the conditional probability 
distribution lE^v-i I Sjvl, which specifies the linear operator Rn-i Rn in (56), ap- 
pears to be a striking reduction of the long-range nature of interactions in the underlying 
Manhattan EDM system.^ 

Remark 3. The "localization" of the original problem above has been achieved through 
its connection with the product moments of the auxiliary GPRF and the possibility of their 
recursive computation (in terms of the conditional product moments) due to the specific 
Markov structure of the random field. A similar line of reasoning, employing the aggre- 
gated level sets and associated conditional product moments, is applicable in three dimen- 
sions for the trivariate GPRFs from [22]. This suggests that such an indirect approach to 
the monomer-dimer problem and its extensions may contain a hidden resource which is 
worth exploring. 

8. One-dimensional Manhattan dimer- monomer system. Consider the Manhattan EDM 
system in the one-dimensional case, where the dimer potential (33) reduces to W{z) := 
a\z\ for z e Z, with a > the attraction force parameter. In this case, the auxiliary GRF, 
which becomes a sequence ^ := (^fc)fegz, is a homogeneous Gaussian Markov chain with 
mean /i from (20) and covariance function cov(^j,^fe) = p'^^*^', where, in accordance 
with (31) and (34), 

p:-e-"^ (60) 
and the shorthand notation (12) is used. In view of (23), the MLE of the sequence takes 
the form 

where 

Af-l 

Mat := Mo:Ar_i(0 = E [] (62) 

k=0 

is the product moment over the discrete interval : — 1, with a : b denoting the set 
Z P|[a, b] of integers from a to b. For convenience of calculations, we will use a centered 
sequence a :— (o'fc)fcgz, defined by ak Cfc ~ P ™d satisfying an autoregressive equation 

cTfc+i = pc^fe + (63) 

for k ^ 0. Here, the initial condition do and the innovations wo,wi,W2, ■ ■ ■ are indepen- 
dent standard normal random variables, and 



7 



Vl-P"- (64) 



*Note that similar recurrence relations, employed in the statistical mechanical transfer matrix method [3] and 
its corner version [2], are known for spin systems with short-range (for example, nearest neighbour or round-a- 
face) interactions. 
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The autoregressive sequence a can be obtained by sampling an Ornstein-Uhlenbeck dif- 
fusion process Y :— (Ft)tgii+ with time step e := a/3 as ak := Ifec- The process Y 
is governed by an Ito stochastic differential equation dYt — —YtAt + v^dWt, where 
W := (yVt)tgR^ is a standard Wiener process, independent of the standard normal initial 
value io- 

9. Conditional product moments and transfer operator. For the one-dimensional Man- 
hattan EDM system of the previous section, the product moment (62) of the auxiliary 
Gaussian Markov chain ^ is representable as 

Mjv = EgAr(ao). (65) 

Here, for any ^ 0, the function Qn : M M is defined by the conditional product 
moment 

QA,(ao) := Mo:Ar_i(C ko) = E ( [] 4 I ao J , (66) 

\ fc=o / 

with Qo = 1- Note that Qat is a univariate counterpart of the function associated 
with the two-dimensional Manhattan EDM system by (54). Thus, by employing a line 
of reasoning similar to Theorem 7.1, and using the Markov property and homogeneity of 
the sequence ^, it follows that the function Qn+i is expressed in terms of Qn from (66) 
through a linear transfer operator G as 

Qw+ilco) = Mo:jv(C I o-o) 

= ^oE(Mi:^(^ I (To, (Ti) I fio) = eoE(Mi^^(C I ai) | ao) 

= (/i + ao)E(QAr(ai) | do) G(Qw)(ao). (67) 

Here, the last expectation is taken over the transition probability law of the sequence a 
which, in view of (63) and (64), is described by 

IfTi I (Tol =AA(pao,7'), (68) 

that is, the conditional Gaussian distribution with mean E(ai | co) = puo, variance 
var((Ti I (Tq) — 7^ and PDF 

Pi|o(y I x) e-(«-^-)'/(27')/(\/2^7)- (69) 
By induction, (67) and (68) imply that the function Qjv is a polynomial of degree N with 
the leading coefficient limj;^oo(QAr(a;)/a;^) = Y\^Zo — The latter also 

employs the property that if X ~ Af{m, s^) with a fixed variance s^, and L is a polynomial, 
then EL{X) is a polynomial in the mean value ni with the same leading term as L. In 
particular, the first three of the polynomials Q n are computed as 

Qlio-o) = + (To, 

Q2(o'o) = (m + o-o)E(Ai + (Ti I (To) = (Ai + o-o)(Ai + po-o), 
Q3(o-o) = (m + o-o)E((m + (Ti)(Ai + pcTi) I (To) 

which shows that the recursive calculation of their coefficients, which are needed for (65), 
in the standard basis of monomials l,x,x'^, . . . quickly becomes complicated. In Sec- 
tion II , we will perform these calculations in a more suitable basis of Hermite polynomials 
after establishing the boundedness of the transfer operator G which governs the recurrence 
equation (67). 
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Remark 4. A similar recursive computation of conditional product moments can also be 
developed for a wider class of homogeneous Gaussian random sequences with rational 
spectral densities which are realizable as a hidden Markov chain with a higher dimensional 
state space. 

10. An upper bound for the moment Lyapunov exponent. Let H denote the Hilbert 
space of functions /i : M M which are square integrable over the standard normal PDF 
V from (8) and are endowed with the norm {h, h) and inner product (/, g) := 

f {x)g{x)v{x)(ix = 'E{f{X)g{X)), where X is a standard normal random variable. 
The Hermite polynomials (9) form an orthogonal basis of the space H, with ||-fffc|| = V^; 
see, for example, [24, pp. 19-20]. The following theorem establishes an upper bound for 
the II • 1 1 -induced norm of the transfer operator G defined by (67): 

|||G|||:= sup ^1^. (70) 

Theorem 10.1. The norm (70) of the transfer operator G in (67), associated with the one- 
dimensional Manhattan EDM system, satisfies 

1 4- n-2. 

(71) 

Proof. Let /i € H. Then, by employing a change of measure technique, it follows that for 
any real x. 




G{h){x) = (fi + x) h{y)piioiy \ x)dy 

= i^i + x) h{y)q{x,y)iy{y)dy={p + x){q{x,-),h). (72) 



Here, 

, , Pi\oiy I poi{x,y) 
iy{y) v{x)v[y) 

with pi|o the transition PDF of the sequence a described by (69); v is the standard normal 
PDF from (8), and 

Poi(x,2/) e(2p-^--'-^')/(27^)/(27r7) (74) 
is the joint PDF of co and ai. Application of the Cauchy-Bunyakovsky-Schwarz inequality 
to (72) yields 

{G{h)(x)f ^{^l + x)^q{x,■)nhf. (15) 
An upper bound for the transfer operator norm (70) can now be obtained by integrating 
both sides of (75) with respect to the standard normal distribution and dividing the result 

hyWhf: 



Gir ^ / l^{x){^Ji + xY l^J^iQix, y))My)dy J dx 

{fi + x)^s{x^ y)dxdy. (76) 

Here, in view of (73), (74) and (64), 

six,y) := {q{x,y)fiy{x):y{y) = ^^"i^^'^)-' 

iy(x)iy[y) 

1 f ipxy - {1 + p^)x^ - {l + p^)y^ \ 
- exp — . (77) 
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The function 7^5 : IR+ is recognizable as a bivariate Gaussian PDF with zero mean 

and covariance matrix 



1 + p2 -2p ■ 
-2p 1 + p2 



1 



2p 1 + p^ 



whose determinant is equal to one and the diagonal entries describe the common marginal 
variance (1 + p^)/7^. It is the latter quantity that enters the right-hand side of (76) through 
(77) as 



IIIGII' ^ \ [ l^s{x, y){fi + x fdxdy = \ ( f,' + 



^2 



thus proving (71). □ 

Now, by using the Cauchy-Bunyakovski-Schwarz inequaUty and the submultiplicativity 
of the operator norm ||| • |||, it follows that the product moment (65) satisfies 

Afjv = {Qn, 1) IIQa.I1 - I1G^(1)|| \lGf (78) 

for any N ^ 0, where 1 denotes the identically unit function. Hence, the MLE (61) admits 
an upper bound m(^) ^ In |||G|||, which, in combination with (71) and (64), yields 

2 l-p^)2 



m(Os^;^ln + -;rln(l-p^). (79) 



I(o-o,o-i) 

Remark 5. The last term in (79) is recognizable as Shannon's mutual information I(cro, <ti) 
between the random variables uo and cti (see [9, 17]), which, in the Markov case being 
considered, coincides with I(tTi,cr^o)> where cr^o '■= {o'k)k^o is the past history of the 
sequence a up until time 0. 

Remark 6. The proof of Theorem 10.1 shows that the inequality (76) remains valid for 
any homogeneous Markov random sequence ^, and hence, its MLE, in this more general 
case, admits an upper bound 



2 Jr2 v(x)v{y) 



where v and poi are understood as the invariant PDF of ^ and the joint PDF of ^0 and ^1, 
respectively. 

1 1 . Matrix representation in the basis of Hermite polynomials. As an orthonormal ba- 
sis of the Hilbert space EI from Section 10, we will use the functions /iq, /ii, /i2, • ■ • obtained 
by scaling the Hermite polynomials (9) as 

hu := Hu/^\. (80) 

Their significance for computing the product moment (62) is that, in view of the leftmost 
equality from (78), 

Mn = qNfl. (81) 

where 

qN.k ■= {Qn, hk) (82) 
are the Fourier coefficients of the polynomial Qm with respect to the orthonormal basis 
(80), so that 

N 

QN^^qN,khk- (83) 

fe=0 
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We will show that the coefficients qn.O: ■ ■ • j Qn.n satisfy a recurrence equation in TV ^ 
which can be "encoded" in an equivalent matrix representation of the transfer operator 
G defined by (67). To this end, we will employ a lemma below which, although it fol- 
lows from the properties of the Ornstein-Uhlenbeck operator [32, Proposition 1.5.4(v) on 
p. 233], is provided with a proof for completeness of exposition. 

Lemma 11.1. Let X be a Gaussian random variable with mean m and variance < 1. 
Then, for any nonnegative integer k, 



EHk{X) = u^Hkim/u), 



(84) 



Proof. By combining the generating function (10) of Hermite polynomials with the mo- 
ment-generating function Ec'^^ = e™^+* ^ of the Gaussian distribution J\f{m, s^), it 
follows that 

+°° ,,k 



fc=0 



^y-(i-s^)r/2 



+00 

E 

fe=o 



{uyf 
kl 



Hk{m/u), 



which holds for all y. Comparison of the leftmost and rightmost sides of this identity yields 
(84). □ 

Now, let R, A and be linear operators on the Hilbert space H which act on the basis 
functions (80) as 

R{hk) := p'^hk, A{hk) := Vkhk-i, A^hk) := Vk + lhk+i. (85) 

Both A and ^4^^ are defined (and are mutually adjoint) on the linear manifold of those 
functions h E M whose Fourier coefficients in the basis (80) are square summable with an 
appropriate weight: J2k>i ^(^^ ^k)^ +oo. The action of the operators A and in (85) 
on the basis functions is identical to that of the annihilation and creation operators on the 
eigenfunctions of the Schrodinger equation for the quantum harmonic oscillator; see, for 
example, [35, pp. 52-53] and [42, p. 91]. In this sense, A + A"^ (which is a symmetric 
operator) corresponds to the quantum-mechanical position operator [24, pp. 210-211]. 
Although A and A^ are unbounded, their compositions AR and A^ R with the operator R, 
as well as R itself, are bounded operators on the whole space M, since < p < 1 in view 
of (60). 

Theorem 11.2. The transfer operator G in (67) can be expressed in terms of the operators 
R, A and A^ from (85), and its matrix representation in the basis (80) is given by 



G = {^iI + A + A'')R = 





p 
















1 


pp 


V2p^ 











* 





V2p 


pp^ 


^p' 
















Vip' 


pp^ 


4 

pp^ 




































pp"^ 






* 


* 




* 


* 





(86) 



Proof. In view of (68), application of Lemma 11.1 with m pa^ and s 7 from (64) 
yields the identity 

E(Hfc(ai) I do) -p'=fffc((To), 

by which the random sequence [p^-'^ Hk{<Jj))j^Q is amartingale with respect to the natural 
filtration of the homogeneous Gaussian Markov chain a for any given fc ^ 0. Hence, the 
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action of the transfer operator G on the fcth Hermite polynomial is described by 

G{Hk){cro) = (Ai + ao)E(i?fe(ai) | ^o) =p''{f^ + cro)Hk{cro) 

= p''i^^Hk{ao) + kHk-li<J„)+Hk+li<Jo)), (87) 

where we have also used the recurrence relation Hk+i{x) = xHk{x) — kHk-i{x). Upon 
division of both parts of (87) by ^/Itl, the result is representable in terms of the basis func- 
tions (80) and the operators (85) as 

G{hk) = p^{phk + ^khk-i + \/k+lhk+i) 

= fiR{hk) + A{R{hk)) + AHRihk)). (88) 

Since the last relation holds for any fc ^ 0, the transfer operator G is indeed expressed in 
terms of R, A and A^ as in (86). The first of the equalities (88) yields the kth column of 
the matrix representation of G in the basis (80) on the right-hand side of (86). □ 

Since < p < 1, the matrix representation (86) implies that G is a compact operator. 
With a/R denoting the square root of the operator R from (85) defined by ^/R{hk) — 
p^/'^hk, the spectrum of G coincides with that of a self-adjoint operator ^/R{fj,I + A + 

) \/R and is all real. Moreover, since p, > 0, the convex cone E1I+ of functions /i e H 
with all nonnegative Fourier coefficients {h, hk) with respect to the orthonormal basis (80) 
is invariant under G. Hence, by the Ruelle-Perron-Frobenius theorem, the operator G has 
a unique (up to a positive scalar factor) eigenfunction 

+00 

Q{x):=J29khk{x) (89) 

fc=0 

in H+, which corresponds to the largest eigenvalue A of G. By applying (86) to the Fourier 
coefficients (82) of the polynomial Q^, it follows that for any N ^ 0, they satisfy the 
recurrence equation 

QN+l.k - Vkp''-^qN,k-l + Pp'^qNM + Vk + lp''+\NM+l, s$ fc < TV + 1. (90) 

Here, qpf^k — for k > N and fc < in view of (83), and the recursion is initialized 
by qo,a = 1 and go,fe = for any k ^ 0. Since the cone H_|_ is invariant under G, a 
straightforward induction yields the nonnegativeness q^ ^ ^ 0. Furthermore, the ratio of 
the free terms (81) of the polynomials Qn and Qn+i converges to the largest eigenvalue 
of the transfer operator G as 

lim *™=A, (91) 

which can be used for finding this eigenvalue numerically as it is done in the power method. 
The dependence of the MLE m(^) = hi A on /i and p, obtained by numerical computation 
of the eigenvalue A, is presented in Fig. 4. 

12. Connection with a pantograph equation eigenvalue problem. We will now discuss 
an alternative representation for the eigenvalue A from (91). Consider an exponential type 
generating function of the Fourier coefficients (82): 

k=o ^ k>a ^K. 




^ ^-Hk{-) ) = e-''/'-EiQNiX)e'''), (92) 
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m(0 




Figure 4. The dependence of the MLE m(^) of the auxihary autore- 
gressive Gaussian random sequence ^ for the one-dimensional Manhat- 
tan EDM system on the parameters and p from (20) and (60). Note that 
m(^) — >■ In /J, as p — )■ 0+, and m(^) — +oo as p — 1—. Another limit, 
lim^_,.o+ ni(^), which is a function of p, is unknown. 

where the second sum consists of finitely many nonzero terms (with k ^ N), and X is a 
standard normal random variable. Here, use has also been made of the generating function 
of Hermite polynomials (10). The definition of Tjv implies that it is a polynomial of degree 
N. 

Theorem 12.1. The polynomials Tpf, defined by (92), satisfy the recurrence equation 

Tjv+i(t) = {p + t)TM{pt)+ pT'^{pt) 

= e-(^+*)'/2at(rAr(pt)c(^+*)'/2) K{TN){t) (93) 

Proof. Although the relation (93) follows from (90), we will derive it in a more instructive 
fashion using a change of measure. By combining the definition (92) with the recurrence 
relation (67) and repeatedly using the tower property of iterated conditional expectations, 
it follows that 

= e-*'/2E((p + ao)c*^''E(Qjv(ai) | ao)) = e~'" '^^{{^jl + (To)e*"''Qiv(^7i)) 

= e-*'/2E(g^(ai)E((Ai + ao)e*"" | ai)). (94) 

Since the Gaussian Markov chain a is reversible, the conditional probability law of do, 
given (7i, is obtained from |(Ti | ctoI by swapping o-q and cti, so that |(To | cil = 
J\f{p(Ji,j'^), with 7 given by (64). Hence, the conditional moment-generating function 
of (To is 

Otiai) := E(e*"« | cti) = e*'"^i+*'^'/2_ (95^ 

Now, since 

E((Ai + (7o)e*"° I ai) = {p + dt)et{ai), (96) 
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then substitution of (95) and (96) into the right-hand side of (94) represents the latter equa- 
tion as 

= e-*'/2(^ + dt)ie'"/^TN{pt)) = + t)TNipt) + pT'j^ipt), (97) 
thus establishing (93). In (97), use has again been made of (64) and (92). □ 

The linear operator K, defined by (93), is the composition of a linear first order differ- 
ential operator with an argument scaling operator. Since K consists in differentiation and 
scaling, it is easier to iterate than the integral operator G in (67) which involves conditional 
expectation. With the initial condition Tq = 1, the subsequent three of the functions Tjv 
are computed as 

ri(i) = ^i + t, 

T2{t) = i^L + t){fi + pt) + p, 

Tsit) ^ itl + t){{^I + pt){p + ph) + p} + pi^I + ph + p{fl + pt)). 

By the relation Tjv+i(0) = pTn{0) + pT'j^{0), which follows from (93), the limit (91) is 
representable as A = lim^^+oo(T;v+i(0)/rAr(0)) = fi + plimjv^+^(r;V(0)/TAr(0)). 
Here, the rightmost ratio, which is the logarithmic derivative (lnrAr(i))'|(=o, is equal to 
<lN,i/qN,o- The generating function of the Fourier coefficients of the eigenfunction (89), 
defined by 

T(i):-^qfc — , (98) 

fe=o 

is a solution of the eigenvalue problem K{T) = AT for the functional-differential operator 
K, that is, (/i + t)T{pt) + pT'{pt) = XT{t). Upon introducing a rescaled independent 
variable r pt, the eigenvalue problem takes the form 

r'(r) = -Tir/p) t^±^T{r), (99) 
P P 

which is a pantograph equation [5, 10, 26, 31, 44] with a variable coefficient. The presence 
of the independent variable r together with its scaling r/p as arguments of the unknown 
function T makes this functional-differential equation hard to solve, let alone finding A as 
the largest eigenvalue for which the pantograph equation (99) has a feasible solution (98). 

13. Concluding remarks. We have established a novel connection between the product 
moments and moment Lyapunov exponents of homogeneous Gaussian random fields on a 
multidimensional integer lattice with the partition function of an elastic dimer-monomer 
system which extends the classical monomer-dimer problem of equilibrium statistical me- 
chanics by allowing for long-range interactions. A salient feature of the auxiliary GRF, 
which encodes the thermodynamic properties of the EDM system, is that its mean value 
and the covariance function are the Boltzmann factors, associated with the energetics of 
the underlying particle system, and, therefore, are both nonnegative. Any insight into how 
the MLE of such a GRF depends on its spectral density function in an arbitrary dimen- 
sion would shed light on the solution of the EDM problem, including its particular case, 
the classical monomer-dimer problem, which remains unsolved in dimensions three and 
higher 

To this end, we have outlined a research paradigm which builds on the observation that 
a specific spatial Markov structure of the GRF, present, for example, in Gaussian Pickard 



MOMENT LYAPUNOV EXPONENTS OF GAUSSIAN RANDOM FIELDS 



23 



random fields, facilitates recursive representation of the product moments of such GRFs. 
We have discussed the possibility to recursively compute the product moments with a local 
transition from one level set to another, for a Gaussian PRF, which is the auxiliary GRF 
for the Manhattan EDM system. This paradoxical "localization" of the effect of long- 
range interactions suggests that the product moment approach, proposed in the present 
study, is a promising resource towards the solution of the monomer-dimer problem and its 
generalizations. 

In addition to its links with the statistical mechanical lattice models of interacting parti- 
cle systems, the problem of computing the product moments and MLEs for homogeneous 
GRFs appears to be of interest in its own right from the probability theoretic point of view. 
In this regard, even the deceptively simple case of univariate Gaussian Markov chains re- 
veals rich algebraic and probabilistic structures. In particular, we have linked the product 
moments and MLEs in this case to the ladder operators of the quantum harmonic oscil- 
lator and the pantograph functional-differential equation by using a change of measure 
technique. 
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